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Abstract 

Using a generalized transition probability, it is shown how nonadiabatic calculations, within the 
Wigner-Heisenberg representation of quantum mechanics, can be reliably extended to far longer 
times than those allowed by a primitive sampling scheme. Tackling the spin-boson model as a 
paradigmatic example, substantial numerical evidence is provided that the statistical error can be 
reduced for a wide range of parameter values. 
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I. INTRODUCTION 



Many quantum processes can be modelled by the time evolution of a relevant subsystem 
coupled to a bath. In certain representations, the dynamics of the subsystem can be un- 
derstood in terms of the transitions between quantum states as a result of the nonadiabatic 
interaction with the bath. Various methods have been proposed to formulate numerical al- 
gorithms for nonadiabatic dynamics UQ. The main theoretical difficulty is associated with 
a rigorous description of the back-reaction [5] of the subsystem on the bath (or a systematic 
way to derive controllable approximations of this phenomenon). 

A systematic approach for deriving approximations of nonadiabatic quantum dynamics 
has been proposed some time ago jsj: One can start from the Heisenberg equations of motion 
of the system and perform a partial Wigner transform [6( over only the degrees of freedom of 
the bath. In such a way, one obtains what may be called Wigner-Heisenberg representation 
of quantum mechanics p|. For harmonic baths, such a Wigner-Heisenberg representation 
leads to an exact formulation of quantum dynamics. For more general baths, a linear 
approximation of the propagator leads to a hybrid quantum-classical representation jg]. 
In this case, time evolution is dictated by what is known as quantum-classical Liouville 
equation [9|. This equation generates a unitary dynamics for all the degrees of freedom 
in the system. Integrating over the phase space degrees of freedom [lo|], one can obtain 
stochastic equations for the remaining degrees of freedom satisfying the quantitative criteria 
analyzed in llj|. One of the appealing features of the Wigner-Heisenberg formulation of 
mixed quantum-classical systems is that it leads to a rigorous formulation of their statistical 



mechanics 



121 ] . In both cases, controllable approximations lead to a consistent formulation of 



nonadiabatic dynamics, which is energy- conserving and correctly describes the back-reaction 
mentioned above. From the computational point of view, this requires one to evaluate 
statistically the propagator in terms of a stochastic process describing the occurrence of 
quantum transitions: Wigner-Heisenberg quantum mechanics can be represented in terms of 
piecewise adiabatic trajectories of the bath coordinates, interspersed by stochastic transitions 
between the energy levels of the subsystem. Despite the elegance of such an approach, the 
growth in time of the amplitude of the statistical fluctuations has so far somewhat limited 
the length of the time interval that can be explored. In order to tame this problem two 
approaches must be combined. The first is to develop more accurate approximations for the 
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short-time expression of the Wigner-Heisenberg propagator [13j. This is critical especially for 
high values of the friction. The other is to improve the time-dependent stochastic sampling 
of nonadiabatic quantum transitions [l^j]. Here we are concerned with this second issue. 

In the following, starting from a simple approximation of the short-time Wigner- 
Heisenberg propagator, wc discuss in detail how a sampling scheme recently introduced 
in can dramatically improve the numerical stability of the calculations at long time (de- 
pending upon the numerical parameters, the time span can be up to an order of magnitude 
longer than in previous schemes). Our sampling method is based upon a suitable generaliza- 
tion of the transition probability that effectively filters some trajectories which would cause 
an excessive energy fluctuation, within an approximate representation of the back-reaction. 
Such a filtering limits the growth in time of the statistical fluctuations. This in turn allows 
the numerical method to access far longer integration times. 

The structure of this paper is the following. In Sec. [Ill we briefly sketch the theory of 
the Wigner-Heisenberg representation of quantum mechanics. Section II III illustrates how 
nonadiabatic dynamics can be derived within the Momentum- Jump (MJ) approximation. 
In Section HVl we introduce a generalized sampling scheme for the nonadiabatic propagator. 
In Section [V] we present the results of some numerical calculations on the nonequilibrium 
dynamics of the spin-boson model. Finally, in Sec. I VI I we discuss our achievements and 
propose future directions. 



II. WIGNER-HEISENBERG REPRESENTATION OF QUANTUM MECHANICS 



Consider a system defined by the total Hamiltonian operator 

H = H s + H B + H S b , 

where the subscripts S, B and SB stand for subsystem, bath and coupling. 
The Heisenberg equation of motion of a generic operator \ can be written as 

H 



(1) 

respectively. 
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where 3 C is the symplectic matrix 
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We assume that the bath Hamiltonian depends upon a pair of canonically conjugate 
operators, X = (R,P), and that the coupling has the form H SB = H S b{R)- A partial 
Wigner transform can be introduced as 

Xw(X) = J dze lP ^ h (R - ^\x\R+ Z ~) , (4) 
for the generic bath-dependent operator \ an d, analogously, for the density matrix p as 



Pw{X) 



(2nh) 



1 J^p-'-M+i) 



(5) 



where X = (R, P) are canonically conjugate classical variables in phase space. The Wigner- 
Heisenberg equations of motion are obtained upon taking the partial Wigner transform of 
Eq. 0: 



H W {X) 
Xw(X,t) 



(6) 



where 



3 2 U k °kj U J 







(7) 



In Eq. ([7]) we have used the symbols a ^ = a /dX k and d \ = d /dX k to denote the 
operators of derivation (with respect to the phase space point coordinates) acting on the 
right or left, depending on the direction of the arrow. The sum over repeated indices must 
be understood in Eq. ((?]) and in the following. The mixed Wigner-Heisenberg form of the 
Hamiltonian operator, Hw, is 



Hw{X) = H s + H W , B {X) + H W! sb(R) ■ 



(8) 



Equation (jH]) provides a mixed Wigner-Heisenberg representation of quantum mechanics, 
where operators also depend upon phase space (c-number) coordinates. Such a represen- 
tation is completely equivalent to the usual Heisenberg representation, but calculations are 
very difficult in general. However, in the case of quadratic bath Hamiltonians 



W,B 



Ehr + 5 



k=l 



(9) 



where (R^, Pk), k — 1, . . . , N, are the coordinates and momenta, respectively, of a system of 
N independent harmonic oscillators with frequencies Uk, and for interaction Hamiltonians 
of the type 

Hw,sb = V b (R)®H' s , (10) 

where Vb{R) is at most a quadratic function of R and H' s acts only in the Hilbert space of 
the subsystem, Eq. (jBJ) can be rewritten using the antisymmetric operator matrix 



lin 





-1 — — d tBkj~3 j 



1 + 4r d hBhi~3 



2 ^ k^kj U j 





(11) 



Actually, it can be shown that for the class of Hamiltonians specified by Eqs. fl9]) and (TTO 



T> -»■ T> 



lin 



(12) 



holds exactly. For more general bath Hamiltonians and couplings, such a substitution 
amounts to performing a quantum-classical approximation js]. What matters here is that 
for the class of systems in which we are interested, Eq. (|T2"|) is exact and provides via Eq. (j5J) 
a Wigner-Heisenberg formulation of quantum mechanics which can be numerically simu- 
lated [12]. 



III. MOMENTUM-JUMP NONADIABATIC DYNAMICS 

In order to devise a numerical algorithm, Eq. must be represented in some basis. The 
adiabatic basis, defined by 



Hw{X) 



P 2 
2M 



+ h w (R) 



(13) 



and the eigenvalue equation, 



h w (R)\a;R) = E a (R)\a;R) 



(14) 



naturally leads to surface-hopping schemes, as is clarified below. In such a basis the 
quantum-classical evolution is 



act 

Xw 



(15) 
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where 



i^aa',/3/3' ~ («<<W + iL aa /) 5 a p8 a >p> + J a a',/38' (16) 
= iC Q aQ ,8 a fi8 a ipi + J aa ',pp' (17) 



is the quantum(-classical) Liouville operator In turn, the Liouville operator is defined 
in terms of the Bohr frequency, 

u a AR) = EaiR) ~ h EAR) , (18) 



the classical-like Liouville operator for the bath degrees of freedom 

M"dR + 2~\ w+ w ) ' dP 



»'^=-S~+if^+^)-4. ( i9 ) 



and the transition operator 

Jaa'BB' 



P ( I AE a p(R)d af3 (R) d\ 



+ p d* fnJ, 9 1a r 2 m 

" m ■ I, 1 + 2 ' ap 1 ^ ' (20) 



where AE a p(R) = E a (R) — Ep(R). The J operator defined by Eq. f[20]) is responsible for 
the nonadiabatic transitions between the energy levels of the quantum subsystem as a result 
of the coupling to the bath. 

It is in general difficult to devise numerical algorithms that implement the exact form 
of the transition operator in Eq. f[2"Uj) . The Momentum- Jump approximation consists in 
approximating the J operator as 

Jaa',BB' ~ J^a',/38' = T^B^a'B' + 7j^/<L/3 , (21) 

where 

P f "l AE aP {R)d aP {R) d] (22) 



T^b = ^ ■ d aB (R) exp 



2 £ . d aP {R) BP 
In this simplified form, the action of J MJ on the momenta P (the back-reaction) when a 
quantum transition occurs can be calculated exactly. 

For simplicity, in the following we consider that all masses M are equal. When an 
a — > (5 transition occurs, the application of the momentum-jump operator to the momenta 
P produces a shift, P — > P', defined by 

P = P + A^P, (23) 



where 



A™P = -(P-d a0 )d at 



+ d a0 sign(P ■ d aP ) J (P ■ 4/j) 2 + MAE a p , (24) 

and d a p is the unit vector associated to the coupling vector in the multidimensional space 
of all the particle coordinates. The expansion of the square root on the right hand side of 
Eq. f[24"]) provides the approximated momentum-shift 

~ MJ l AE a p(R) j 

A *0 F = » P o — • (25) 

The exact momentum-shift, given by Eq. (1241) . conserves the energy in eve ry q uantum tran- 
sition. Accordingly there are no problems with so-called "frustrated hops" [l8j] . Instead, the 
property of energy conservation is broken only by the approximated form of the momentum- 
shift given by Eq. ( 125|) . In this work, the numerical propagation of the trajectory is performed 
using the exact energy- conserving momentum- shift rule. In the next section, we will discuss 
the sampling probability for the quantum transitions. Such a transition probability is not 
uniquely fixed. Hence, we can modify it using a filtering scheme based on the approximated 
momentum-jump. We show that such a filtering leads to a substantial reduction of the 
statistical error in the calculations. 

IV. STOCHASTIC SAMPLING OF NONADIABATIC DYNAMICS 

The formulas of the previous section can be implemented through an elegant and simple 
algorithm which is based on a sequential time step expansion of the Dyson propagator [19] . 
For a small time step At the quantum-classical propagator is approximated as 

(e^ MJ ) - e^> (<W<W + AtJ^) . (26) 

V / act', 80' 



( |26|) reproduces 



19( | . However, from 



One can prove that the concatenation of short time steps according to 
exactly the Dyson integral expansion of the operator exp (iAtjC) aa , on, 
a computational point of view, it is impossible to evaluate exactly the concatenation of the 
short time steps. What one can do is to consider the sequential short time propagator in 
Eq. (1261) as a stochastic operator and sample the action of the transition operator J according 
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to suitable transition probability. Within such a stochastic evaluation of the Dyson expan- 
sion of the propagator, the action of the operator iC^^p, is that of generating a piecewise 
adiabatic process for the momentum and coordinate variables, P and R, respectively, over 
energy surfaces identified by the state labels a and a'. 

The practical implementation of the stochastic propagator requires the sampling of the 
transitions. A primitive choice for the transition probability is given by: 

^ (x ' At) -i + liH^)|At' (27) 

which also determines the probability of not making any transition in the time interval At 

as 



1 



l + \£-d«p(R)\At ' (28) 

One can exploit a certain arbitrariness in the definition of the transition probabilities, 
Eqs. ()27p and f[28]) . in order to define a filtering of the nonadiabatic transitions. As we 
illustrate by some numerical calculations, this reduces the statistical noise at long time. 
To this end, we consider the variation of the energy, as calculated from the approximated 
momentum shift in Eq. (j25p . After the quantum transition one has 

«* = Q + em - g + (29) 

where 

P' = P + AP afj (30) 

are the new momenta arising from the quantum transition a — > fl. We define new transition 
probabilities as 

T> (Y M\ At\^-d a p{R)\u{c £ ,£ a p) 
l + At\-jg ■ d a f3{R)\uj{c £ ,£ al3 ) 

and 

Q a p(X, At) = l-Vaf, 

P 3 Wtt: , ^ > ( 32 ) 



1 + At\^ ■ d af s(R)\u(c £ ,£ al3 ) 
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where 

{1 if £ a p < cs 
(33) 
otherwise 

The numerical parameter eg controls the amplitude of the energy fluctuations in a transition. 

Eq. (l3Tj) defines a generalization of the primitive sampling scheme. The analytical form of 
the weight a; is to a large extent arbitrary. The choice adopted in the calculations reported 
in this paper and in the previous jl4] is devised by requiring that the energy fluctuations 
caused by the approximated momentum-shift rule (1251) are not too large. In what follows, 
we illustrate by means of numerical calculations that the particular choice given by Eq. ( 133]) 
leads to a dramatic reduction of the statistical noise of the algorithm. 

V. NUMERICAL EXAMPLE 

In order to illustrate the effectiveness of our generalized sampling scheme, defined by 
Eqs. (I31ti33p . we report the results of calculations performed on the nonequilibrium quan- 
tum dynamics of the spin-boson model. Such a model is defined (using adimensional coor- 
dinates [17j) by 

H s = ~na x (34) 

N 

H w ,sb = -Zz^CjRj , (35) 

3=1 

where a x and a z are the spin Pauli matrices and Cj are numerical coefficients whose value 
is fixed by requiring that the spectral density of the system has a Ohmic form. The bath 
Hamiltonian is given by Eq. fl9]). For this model the adiabatic basis is known analytically j^oj]. 
One can think of situations in which at t = the spin is excited in state up and the quantum 
harmonic modes are at thermal equilibrium as if there were no coupling before t — 0. For 
later times, the coupling can be switched on and one can calculate, for example, the d ecay 



of the population, (a z (t)). The dynamics of the spin-boson model has been well studied 20], 
and thus provides us with a good system to check the efficacy of our generalized sampling. 

Here we compare the results of calculations performed with the generalized sampling 
scheme to those of the primitive algorithm. In all calculations a maximum number of 
n = 2 nonadiabatic transitions per trajectory has been considered. For each matrix element 
propagated in time, the phase space ensemble has been composed of a number N mcs = 



2.5 x 10 5 points. The integration time step of the phase space trajectory has been taken 
as dt = 0.01 in dimensionless units [17]. In Figs, [fl |2] and [3] we display the results for 
the observable (a z ) against time for three sets of parameters. Figure Q] shows the results 
obtained for the parameters ft = 3, Q = 1/3 and £ = 0.1. The generalized sampling with 
a parameter c e = 0.05 stabilizes the dynamics over a time-interval four to five times longer 
than that achieved with the primitive sampling scheme. Figure [2] displays the results for 
(3 = 0.25, Q — 1.2 and £ = 2, which represents higher friction for the spin-boson dynamics. 
In this case, using c e = 0.1, the generalized sampling allows one to extend the dynamics for 
at least a twice as long a time interval than that of the primitive sampling. Longer sampling 
is also achieved for /3 = 1, Q = 0.4 and £ = 0.13 with c e = 0.1. Results are displayed in 
Fig. [3] 

The introduction of the weight w, defined in Eq. (13"3"|) . in the generalized transition prob- 
ability given by Eq. (I3"T|) realizes a pruning in the ensemble of the sampled quantum tran- 
sitions, while keeping the important part of the nonadiabatic effect. Figure H] shows the 
fraction of trajectories in the ensemble that performed n — 0, 1, 2 nonadiabatic transitions 
against time, for the generalized sampling (main figure), and the primitive sampling (inset). 
Numerical fluctuations are reduced since the generalized weight makes Q = 1 when a transi- 
tion is rejected. Our calculations for the nonequilibrium dynamics of the spin-boson model 
show that the pruning allows one to sample longer times in a stable way, while reproducing 
the short-time numerical results of the primitive algorithm, provided a suitable choice of c e 
is made. 

These results, together with the calculations at lower friction reported in jlj], provide 
substantial numerical evidence that the generalized sampling is able to reduce the statistical 
error at long time. 
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FIG. 1: Comparison of primitive (A) and generalized (•) sampling for f3 = 3, = 1/3, and £ = 0.1. 
A parameter value of eg = 0.05 was used. Two quantum transitions per trajectory were included. 
The error bars associated the primitive sampling start growing from t = 10 and become enormous 
after such threshold. The error bars associated to the generalized sampling become distinguishable 
from the bullet only after t = 16 and remain, nevertheless, very small up to t = 20. 

VI. DISCUSSION AND CONCLUSION 



We have shown how our generalized sampling scheme can improve the numerical stability 
of nonadiabatic dynamics at long times. For devising our sampling scheme, we have been 
guided by the principle of energy conservation. In this paper, our choice of the weight w 
reflects such a line of thought. The results of our calculations illustrate in many relevant 
cases how big the reduction of the statistical error is at long time. However, other possible 
choices of the weight w (for example, better suited or tailored to the problem at hand) are 
open for investigation. The important idea is trying to get a sequential weight equal to 
one (most of the times) when an adiabatic segment of trajectory is propagated. This limits 
the growth in time of the statistical fluctuations induced by the sampling of nonadiabatic 
dynamics. One can speculate that other choices of w can be tailored and optimized for 
the specific process investigated. We stress that these alternative choices of w will all be 
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FIG. 2: Comparison of primitive (A) and generalized (•) sampling for (3 = 0.25, O = 1.2, and 
£ = 2. A parameter value of Cf = 0.1 was used. Two quantum transitions per trajectory were 
included. The picture shows clearly the error bars associated with the primitive sampling. Those 
associated with the generalized sampling are pretty much indistinguishable from the bullet symbol 
up to t = 8. 

encompassed by the general structure of Eq. f[3"Tl) . 

We must remark that in order to reduce the systematic error at high friction values, one 
is forced to use a higher-order approximation of the short-time propagator, which has been 
first introduced in [13 |. Here, our concern has been only the reduction of the statistical 
error. To this end our sampling scheme is very successful. As a matter of fact, the accurate 
calculations of [lsj have been obtained with a simple filtering recipe, which basically amounts 
to a redefinition of those values of the observable which are too high as equal to a threshold 



constant 



21] . We believe that the systematic application of our computational theory can 



start an exciting time when the power of nonadiabatic algorithms, stemming from the work 
of 8[], can be reassessed and further extended by combining our enhanced sampling scheme 
with higher order approximations of the short-time propagator 13]. 

It is also worth mentioning that our sampling scheme can be applied to a more advanced 
version of Wigner-Heisenberg quantum(-classical) mechanics (which has been first proposed 
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9 10 11 12 13 14 15 16 17 18 

FIG. 3: Comparison of primitive (A) and generalized (•) sampling for (3 = 1, Q = 0.4, and 
£ = 0.13. A parameter value of eg = 0.1 was used. Two quantum transitions per trajectory were 
included. The inset shows the short-time behaviour where the two sampling schemes provides very 
close numerical results. The main figure shows the long-time behavior. The error bars associated 
to the primitive sampling are clearly seen. Those associated with the generalized sampling are 
barely distinguishable from the bullet symbol only after t = 16. 



m 



22j and further developed in 23J) that can address time correlation functions with quan- 



tum corrections to the structure of the bath coordinates even with anharmonic potentials. 
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